#!/usr/bin/env python
"""
Fig18A.py

Reproduces figure 18 A in Jonas et al., 1993
"""

import numpy as np
import matplotlib.pyplot as plt

from CA3.library import jonas1993
from neuron import h

cell = jonas1993()

# initial parameters
# voltage-clamp at the soma
VC_patch = h.SEClamp(0.5, sec = cell.soma)
VC_patch.rs = 0.1
VC_patch.amp1 = -70

# set one synapse originally at the soma
# for parameters check page 624 and 649
syn = h.synapse(0.5, sec = cell.soma)
syn.tonset = 0 
syn.tau0 = 0.2 # in ms (tau onset in Jonas et al., 1993)
syn.tau1 = 2.5 # in ms (tau decay in Jonas et al., 1993)
syn.gmax = 300e-6 # 300 pS in Jonas et al., 1993


def simulation(tstop, with_time = False):
    """
    runs the simulation and returns the current and
    time vectors as Numpy arrays
    """
    h.load_file('stdrun.hoc')
    h.v_init = -70
    
    h.tstop = tstop
    VC_patch.dur1 = tstop

    # define vectors
    current = h.Vector()
    current.record(VC_patch._ref_i)
   
    if with_time is True:
        time = h.Vector()
        time.record(h._ref_t)

    h.run()

    if with_time is True:
        return (time, np.array(current)*1000.)
    else:
        return np.array(current)*1000. 


# 1. EPSC at the soma, Rs = 0.1 MOhms
time, soma_lowRs = simulation(tstop = 10, with_time = True)

# 2. EPSC at the dendrite, Rs = 0.1 MOhms
# change compartment 
# (move syn to cell.ds[]), which is at a distance of 120 um from soma
syn.loc(0.5, sec = cell.ds[8])
dend_lowRs = simulation(tstop = 10)

# 3. EPSC at the dendrite, Rs = 10 MOhms
# change rs
VC_patch.rs = 10
dend_highRs = simulation(tstop = 10 )

# 3. EPSC at the soma, Rs = 10 MOhms
# go back to soma
syn.loc(0.5, sec = cell.soma)
soma_highRs = simulation(tstop = 10 )


# Plot everything
plt.plot(time, soma_lowRs, 'k')
plt.plot(time, soma_highRs, 'k.')
plt.plot(time, dend_highRs, 'r.')
plt.plot(time, dend_lowRs, 'r.')

plt.xlabel('Time (ms)')
plt.ylabel('Current (pA)')

plt.show()
